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Abstract 

Generalized power asymptotic expansions of solutions to differential equations 
that depend on parameters are investigated. The changing nature of these expan- 
sions as the parameters of the model cross critical values is discussed. An algorithm 
to identify these critical values and generate the generalized power series for distinct 
families of solutions is presented, and as an application the singular behavior of a 
cosmological model with a nonlinear dissipative fluid is obtained. This algorithm 
has been implemented in the computer algebra system Maple. 
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1 Introduction 



Quite frequently physical models depend on parameters and the predictions 
of these models, when confronted with observations, yield constraints that 
fix the parameters or disqualify the model. Hence, for models described in 
terms of a set of differential equations, it becomes of prime importance to 
investigate the dependence of solutions on the parameters. Now, when exact 
solutions to these equations are not available, or even when they are available 
but they have a complex implicit or parametric form, it may become useful 
to obtain asymptotic expansions of these solutions. And when solutions do 
not admit power series expansions with integer or rational exponents (Pusieux 
series), we are led to consider series with more general terms like "generalized" 
powers with real exponents [1] or exp-log terms [2]. Though the algebraic 
issues of generalized power series have already been investigated (cf. [3]), for 
applications the main problem lies in the dependence of the exponents on the 
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parameters, as the ordering of terms and even the nature of the asymptotic 
expansion may change as these parameters reach some critical values. 

In this paper, we will discuss the problems involved in the determination of the 
critical values of the parameters and an algorithm to calculate the coefficients 
and exponents of the generalized power series in the regions of the parameter 
space where real solutions admit such an expansion. 

Several critical steps in the search for solutions in the form of generalized 
power expansions involve calculations with large number of terms and this 
number grows very fast with the size of the ODE, making hand calculation 
inconvenient. On the other hand, as many of the steps in these calculations 
have a systematic nature, use of computer algebra systems appears ideally 
suited. For this purpose, we have developed a set of routines in Maple (for a 
brief review of a previous implementation, see Ref. [4]). 

The plan of this paper is as follows. In section 2, we present the basics of 
generalized power expansions of solutions to a class of nonlinear ODES. In 
section 3, we review the problem of branching of these solutions as parameters 
cross critical values. The algorithm used to find these solutions is presented in 
section 4 and it is applied to an example equation arising from a cosmological 
model in section 5. Finally, the conclusions are stated in section 6. 



2 The generalized power expansion 



We will consider nonlinear ordinary differential equations of the form 



where the coefficients A { and the exponent may depend on parameters 
Pi, . . . ,pq. In the cosmological setting, y is frequently a monotonic function of 
either the scale factor or the Hubble rate in an expanding universe (cf. [5] [6] 
[7] [8]). Hence y(t) > and equation (1), if not algebraic, is still well defined. 

In case that the general solution to (1) is not available, we are usually inter- 
ested in obtaining some information about it in the form of an asymptotic 
expansion for the limits t — > + and t — > oo, of the form 



y(t) 



(2) 
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In homogeneous cosmological models, t is the universal time and these limits 
frequently correspond to the behavior of the solution near the initial singu- 
larity (the "Big Bang") or at large time. Here, Cj and rij are real, in principle 
functions of Pi,...,pq, Cj 7^ and the exponents rij form an ordered set: 
ni < n 2 < ■ ■ ■ for t — > + and ri\ > n 2 > • • • for t — > oo. So C\t ni is the leading 
term, ri\ 7^ 0, t nj+1 /t nj — > in either limit and the set {t nj } constitutes an 
asymptotic scale. Inserting the series (2) into equation (1), and performing 
the necessary asymptotic expansions, we get another generalized power series 



where Ck and are also real functions of pi, . . . ,pq, and the exponents e k 
form an ordered set: t\ < e 2 < • ■ • for t — > + and e\ > e 2 > ■ ■ ■ for t — > oo. 
For the equation (1) to admit a solution with expansion (2), each of the Ck 
must vanish and this set of equations fixes in principle the Cj,rij in (2) up 
to r — 1 of them that remain free and arise from the integration constants of 
the general solution of (1) (the arbitrary constant corresponding to the time 
translational symmetry is fixed to 0). Each set of solutions {cj,rij} yields a 
series representation of a family of solutions to (1). Further, the constraints 
that these Cj ; Tlj clX6 real and the rij are ordered may delimit regions in param- 
eter space where generalized power law solutions are feasible. If these regions 
do not contain the values of the parameters that make physical sense, such 
solutions have to be discarded even when they are mathematically correct. 



3 Case branching by parameter variation 

For simplicity, let us consider that equation (1) depends on a single parameter 
p, so that rij = n,j(p) and Cj = Cj{p). As the exponents {rij} have to be ordered, 
a critical value po of the parameter arises when a pair of consecutive exponents 
become equal, n J+ i(p ) = "^j(po) say Moving the parameter p across po may 
cause changes in the nature of the solution, hence in its series expansion, 
e.g., the asymptotic scale involved, and these effects show up on the behavior 
of the function Vj(p) = rij + i(p) — rij{p) in a neighbourhood of po- Let us 
consider first that Vj{p) is analytic at po so that, excluding very special cases, 
v 5 = dj(p — Po) + • • • holds with some constant dj ^ 0. Hence, the terms 
labeled j and j + 1 on one side of p switch order on the other side. This 
effect is of particular importance, as regards to the analysis of the families of 
solutions, when the leading behavior changes. For the series at po, let us make 
the expansion of this pair of terms 



oo oo 



(3) 



j=l k=l 



Cjt nj + c j+l t ni+1 = t n] [cj + Cj+i + Cj+idj (p — po) hit H ] 



(4) 
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If limcj + i(p — po) 7^ for p — > p 0) as when c J+ i is an arbitrary constant, the 
series picks a logarithmic term, and this case has to be dealt with separately. 
Otherwise, both terms merge into one and a relabeling of terms occur. 

On the other hand, when p is a branching point of Vj as in Vj = dj(p — p) s ^ r + 
• • •, with s odd and r even, an expansion of the real solution as a generalized 
power series exists only on one side and the nature of the solution changes 
when po is crossed. As an example, let us consider the equation 

y + yy + [3y 3 = (5) 



with parameter @, that arises in several cosmological models (e.g., [9] [10] [11] [13]), 
as well as in the analysis of the Painleve equations [12], and in form invariant 
equations [15] [16]. It may be shown that the asymptotic expansion of the 
general solution to (5) is 

y(t) - 7 E c nl n t nr (6) 

1 n=0 



where a± — [1 ± (1 — 8 / 3) 1 / 2 ]/(2 / 9), r = 4 — a is the Kowalevski exponent 
[17], r > for t — > + and r < for t — > oo, c = 1, c n = c n (/3) and 7 
is an arbitrary integration constant. Here, the critical value is (3q = 1/8 and 
Vj — r± — =f8[2(/9o— /9)] 1//2 + - • •. Hence, real solutions admit expansion (6) only 
for (3 < 1/8. For (3 = 1/8, the general solution has a logarithmic expansion 
and for (5 > 1/8, real solutions become oscillatory. 

As the coefficients Cj have to be real and nonvanishing, a critical value p of 
the parameter arises when it is a root of a coefficient Cj(po) = or it is a 
branching point as in Cj(p) = Cj + dj(p — p) s ^ r + ■ ■ •, with Cjo and dj some 
constants, s odd and r even. The expansion (6) to the solution of equation 
(5) shows the second effect as its coefficients are proportional to a, hence they 
are complex for f3 > 1/8. 



4 The algorithm 



We review briefly the algorithm stated in Ref. [4]. The objective of the cal- 
culations is to obtain a truncation of (2) to a finite number of terms, say 
M 

M 

1/m(0 = E^ bj (7) 
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The method of calculation is iterative, so that constants cm,^m, M > 1, 
when not free, are determined by the constants c±, ni, . . . , cm-i, ^m-i that 
were found in the previous steps. 

We start with M — 1, by inserting y 1 = c\t nx into (1). After collecting all the 
terms with the same generalized power, we get a sum of the form 

£>b/i = (8) 

1=1 



We note that, in general, the set of the exponents F = {/z} 1<Ki? is not ordered. 
The i-th term of (1) contributes to (3) with a term of exponent 

9i = it (m - h ) B " ( 9 ) 

h=0 



and coefficient 

E t = AiC^ B * f[(n 1 -h+ 1)^=^ B i = Vi <% (10) 
h=i 



As more than one term of (1) may yield the same exponent (balance), we have 
R < N. The exponents g iy hence the exponents depend linearly on n l5 and 
the pair of terms I and m of (8), to which respectively contribute the terms 
and i'(m) of (1), balance for the exponent n 1 ^ 



^h=o h 


r>h 

B i{l) - 


TDh 

n i'(m)_ 


^2h=o 






TDh 

n i'(m)_ 





(11) 



We note that both the exponents g« and the coefficients Ei may also depend 
on the parameters pi, ■ ■ ■ ,Pq through the exponents Bf and the coefficients 
Ai. In addition, we note that the freedom of multiplying the equation (1) by 
integer powers of y and its derivatives (up to missing solutions where they 
vanish), that change its exponents by B\ — > B\ + m h for some integers m h , 
implies that the exponents in (9) and (10) are representatives of a class of 
exponents related by the transformation laws gi — > gi + YX=o ( n i — ^) mh an< l 
Hi + YJ h =Q mh - 

If the equation (1) admits a solution of the form (2), and yi(t) is its leading 
term, it holds asymptotically y(t) ~ yi(t). Hence on the one hand 

D [y(t)\ ~ dt ei ~ D h t f ^ (12) 
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for some index l\ that we identify by sorting the exponents in the set F. For 
simplicity, only the behavior for t — > + will be considered in the following. In 
this case we get e 1 = f h = min(/ / )i<K R . 

On the other hand we require c\ 7^ and D\ x = 0. As is the minimum of 
the {gi\x <i<N i and we can choose fa 7^ 0, at least a pair of terms in (1) must 
yield this same exponent. Let us say that this minimum occurs for i e /, so 
that C\ = D tl = Ei and this coefficient is a function of C\ and ri\. 

Then, by solving the set of equations fi{n\) = f m {ni), 1 < I < m < R, for 
rii, we obtain the set N e of equality exponents n\ e that make the exponents fi 
exchange order. They delimit intervals within which the sorting of F has to be 
carried out separately. Two cases arise: (a) n\ lays inside any of these intervals, 
(b) rii is an equality exponent. For each interval and equality exponent F 
is sorted and the coefficient of the term with the minimum exponent is 
identified. In case (b) it is verified whether Z\(ci) = is satisfied with a 
nonvanishing root c±, while in case (a) the equation = could determine 
n\ provided that a nonvanishing c\ is feasible. Each pair of real numbers 
(ci, rii) 7^ (0, 0) obtained from this analysis corresponds to the leading term of 
the expansion of a family of solutions. Thus, subsequent calculations to obtain 
higher order terms must proceed separately for each pair. A constant not fixed 
by this procedure corresponds in principle to an integration constant of y(t). 

Additional branching occurs due to the dependence on the parameters. As 
shown in (11), the equality exponents depend on p 1 , . . . ,pg through the B®. 
If a pair of equality exponents exist, n^} and n^} say, and they themselves 
become equal, the equation . . . ,pq) = nf](pi, . . . ,pq) defines a hyper- 

surface in the parameter space (as the exponents Bf may depend on a subset of 
the parameters, only a parameter subspace might need to be considered, and 
when they depend just on a single parameter, hypersurfaces become its critical 
values). Besides, another set of hypersurfaces in the parameter space might 
exist where the equality exponents diverge. Hence, the analysis described be- 
fore has to be done separately inside each of the regions of the parameter space 
delimited by these hypersurfaces, and at their boundaries. 

(h)B li 

For M > 2, when inserted (7) into (1), any factor y M i , h — 0, . . . , r, with 
a nonnegative integer exponent B\ needs to be expanded, or else expanded 
asymptotically up to order M, producing a term of the form Kit^ ni ~ h ^ B (1 + 

\-K M t nM ~ ni ), with some constants Kj. Then, after such expansions, crossed 

terms in the products generate additional terms with larger exponents. For 
example, when M = 2, a term generated by all factors from the leading term 
except one has an exponent g\ — gi + ri2 — n\ > gi. Thus we see that only the 
leading term of (2) can contribute to the leading term of (3). 
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Once the solution coefficients ci, . . . , cm-i and exponents ri\, . . . , um-i (for a 
family of solutions) are determined, the main tasks at step M are: 

(i) Insert yu in (1) and expand (asymptotically to order M) the powers. 

(ii) Collect all the terms with the same power of t. 
Thus we arrive at an expression of the form 

Rm 

D[y M (t)] = J2 D i tfl (13) 
1=1 



where D[ and /; are real functions of pi, ... ,pq, and Rm < NM r+1 as some 
terms in (1) may balance. In general, the sequence of exponents f\, / 2 , . . . , fn M 
is not ordered. As cm and hm only appear in Ck for k > M, the first M terms 
of (13), once sorted after the order of the exponents, are equal to the first M 
terms of expansion (3). In particular, the first M — 1 terms are those already 
found in step M — 1 of the iteration. Then, to identify the candidates for the 
exponent tu we 

(iii) sort the set of exponents F M = {/i} KKflM , and pick the M-th exponent 

flM ■ 

If the equation (1) admits a solution of the form (2), and yuit) is its M-term 
truncation, we get tu = fi M ■ The sorting operation is the most involved part of 
the whole calculation because of case branching. As f t = fi(n M ), those solution 
exponents {n Me } that make a pair of exponents equal, /«( n Me) = /mte) say, 
and satisfy um-i < nMe, separate intervals where a given sorting holds. For 
M > 1 these equality exponents nue arise as solutions to equations of the 
form 

M N r 

E («f - <) % + E E (ft* - 1%) B i = o (1 4 ) 

3=1 i=l h=0 



where the a\ are linear functions of the B\ (see [4] for a geometric interpre- 
tation of this equation). Furthermore, hypersurfaces in the parameter space 
arise as pairs of these equality solution exponents riMe become equal. In (14) 
they enter through the as well as through ni, . . . ,um-i- Hence, step (iii) 
further divides into: 

(iiia) Find the set of the the equality exponents N e = {riMe}- 

(iiib) Identify the hypersurfaces in the parameter space where the equality 
exponents become equal or diverge. 
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(iiic) Sort the {//} for each distinct case. 
The next steps are: 

(iv) Find um (if possible) and Cm- 

(v) Solve Cm = for either cm or um (if not determined in step (iv)). 

A set of routines to deal with steps (iii) and (iv) have been developed in Maple. 



5 Example 



We will show in this section the application of the algorithm sketched in section 
4 to an equation that depends on several parameters and shows some of the 
issues discussed in the previous sections. 

In order to treat dissipative processes in cosmology which are not close to 
equilibrium, a nonlinear phenomenological generalization of the Israel-Stewart 
theory was developed recently [18]. Scenarios in which this kind of processes 
may have occurred include inflation driven by a viscous stress [18,19], and 
the reheating era at the end of inflation [13,14]. In a spatially flat Friedmann- 
Lemaitre- Robertson- Walker universe, Einstein's equations together with state 
and transport equations of the fluid give the evolution equation for the Hubble 
rate H [18] 



k 2 
v 



+ 



3jv 2 
~2a 



' 2k 2 \ E_ 

3 7U 2 ) JJ2 

( ak 2X 



H + 3HH + 



1 - 2 7 N 

, 7 , 



H q - 



7?r 



H 4 [ 
9 



H 2 -* (2H + 3 7 # 2 ) - - 1V 2 H 3 = 



(15) 



where 7, a, v, q and k are parameters describing the thermodynamical prop- 
erties of the fluid. We consider an ordinary viscous fluid so that 1 < 7 < 2, 
0<i><l,a:>0 and k > 0. In the following, we will calculate two term 
truncations of the generalized power expansions of the solutions to equation 
(15) in the limit t — > + , corresponding to the behavior of the solutions near 
the initial singularity. 
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q < 1 


ni < — 1 (7 — g>)ni 
2 

— 1 < n\ < oni — nig — 1 

~2 q ~ 3 

n\ > 3ni — 3 

<7-3 


K q < 3 


ni < — 1 6ni 
ni > -1 3ni - 3 


q > 3 


n\ < — 1 6ni 
2 

— 1 < m < 3ni — 3 

~2 q ~ 3 

n\ > 6ni — nig — 1 

q-3 



Table 1 

Table of leading exponents 
5.1 Leading term 

We begin by inserting the leading term H = c\t ni into (15), looking for solu- 
tions with n\ < and c\ > 0. We get the set of exponents 

Fi = {6 rii, Arii — 2, 5 ni — 1, 3ni — 3, —n\ (q — 7), 6ni — n\ q — 1} 
and the set of equality exponents 

i, — , — , — , --\ 

q — 3 q — 2 g — 4 gj 

These, in turn, are function of q, and we find that there is a critical value 
q — 1 that make them equal. In effect, this is a distinct value as all terms of 
(1) balance and H — c\jt is an exact solution with ci given by the real roots 
of the cubic equation 




2~ V 



97 
T 



7f 




Cl 3 + 



7t! 



Q! 



2A; 2 



Cl 2 + - 1 - 3 - Cl + = 

7 V v^/ 37^f 



This shows that g = 1 delimits different behaviors of the solutions (cf. [20]). 
Besides, as an equality exponent diverges for q = 0,2, 3, 4, sorting of iq must 
be done separately at these values of q and within the intervals they delimit, 
namely (-oo, 0), (0, 1), (1, 2), (2, 3), (3, 4), and (4, oo). The Table 1 shows the 
leading exponent for each case. 
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Thus, we find that the leading exponent switch at n\ = — 1, as well as n\ = 
2/(q — 3) for q < 1 or q > 3, where terms balance. Let us start with these 
equality exponents. For n\ — — 1 and q < 1, we have C\ = D(q — 7) (we denote 
Di by D{fi)) where 



D(q - 7) 



+ 2 



so that c\ = 2/(37) is a leading coefficient. For q > 1, we have Ci = _D( — 6) 
where 



«-0-2(*-5)(H^ 



— V 



- 1 



_2P_ 

37 2 f 2 



so that we get three leading coefficients 
2k 2 2 

Cl = 



37 (k 2 - v 2 ) ' 3 7 (l± v / 2^ 



(16) 



provided that k > v in the first case and v < 1/ \/2 in the third one. For 
the other equality exponent rii = 2/{q — 3) and g < 1 or q > 3, we have 
d = D(3(5 - q)/(q - 3)) where 



8k 2 



q — 3/ a (q — 3) 3 7 t> 2 (q — 3) J 



9 - 1 - - I 



7, 



so that 



Cl 



9 7 3 v A (9-Qq + q 2 ) 
Ak 2 a (7 - 7 g + 2) 



1 

9-3 



(17) 



is the leading coefficient of another family of solutions provided that (9 — 6 q + 
<? 2 )/(7 ~~ 11 + 2) > 0. The remaining cases in Table 1 yield 



D((7-q) ni ) 



9c! 7 - q <y 2 v 2 
2a 



D{Qrii — ri\q — 1) = 3 



Cl 6 i^v 2 ni 
a 
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D(3n! - 

We see that no solution exists in the first two cases, the third case provides a 
solution with n\ = —7/(7 — 1), if 7 7^ 1, and c\ arbitrary, while the fourth case 
shows that D(6n\) — if &; = v or v — l/y/2. For these values, the calculation 
of the leading term has to be done again. 



5.2 Subleading term 



Inserting H(t) = 2/(3'yt) + C2t™ 2 into (15) and expanding the terms with 
noninteger exponents we get the sets of exponents 



F 2 = {-6, 4n 2 - 2, q + n 2 - 6, q + 2n 2 - 5, 5n 2 - 1, 3n 2 - 3, 6n 2 , n 2 - 5, 2n 2 - 4} 



and the set of equality exponents larger than —1, sorted for q < 1 



*SH-5-§.-* 



Thus we find that e 2 = q + n 2 — 6 for n 2 < —q, while e 2 = —6 for n 2 > —q 
and q < 1. In the balancing case n 2 = —q, we find C 2 = D(—6) where 



D(-6) 



4v 2 



37 2 a 



« + 7(9- 2) (y) 9 c 2 



so that a solution exists with the subleading coefficient 
/ 2 V 



a 



c 2 



7(2-9) V3 7/ 



(18) 



For the rest of the cases, we have 



n , fl , v 32 3«2-" 7 - 4+ « (n 2 + 2) V 2 c 2 
^W-6 + n 2 ) = — 



a 
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and neither of them provides a solution. 



Following similar steps, inserting H(t) = c\jt + c 2 t™ 2 , with c\ given by (16), 
into (15) yields e 2 = n 2 — 5 for n 2 < q — 2, while e 2 = q — 7 for n 2 > q — 2. 
At the equality exponent n 2 = q — 2, we obtain the subleading coefficients 
corresponding to the three cases of (16) 



c 2 



Ik 6 - 2 * (k 2 - v 2 ) q ~ 3 v 4 
9^aq (2k 2 - v 2 ) 



8V2v A (V2v =F l) (3 7 /2) 9 (l ± V2i 

faqa {2^2 [± (g — 1) 7 =F 2] f 4 + 2 [7 {2k 2 - l) (q - 1) + 4 (l - A; 2 )" 
[zp 7 (l + 2A; 2 ) (g - 1) ± 2 (4A; 2 - l)] v 2 

(l - 2A; 2 ) (q - 1) - 4A; 2 ] u ± v 7 ^ (g - 1) A; 2 7 }} (19) 



+ 
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Inserting if(t) = cit 2/(9 ~ 3) + c 2 t" 2 , with Ci given by (17), into (15) yields 
e 2 = (13 - 3q)/(q - 3) + n 2 for n 2 < (1 + q)/{q - 3), while e 2 = 2(7 - g)/(g - 3) 
for n 2 > (1 + g)/(g — 3). At the equality exponent n 2 = (1 + q)/(q — 3), we 
obtain the subleading coefficient 



c 2 



3(((g-l) 7 -2)t; 2 + 4fc 2 )7(g-3) 
Ak 2 (g-2)((g-l) 7 -4) 



9 7 3 1; 4 (9 -6g + g 2 ) 
4/c 2 a (7 - 7 g + 2) 



2 

g-3 



(20) 



6 Conclusions 



We have shown some problems that occur when dealing with generalized power 
expansions of solutions to nonlinear ordinary differential equations that de- 
pend on parameters, and we have sketched an algorithm that allows identify- 
ing critical values of the parameters and obtaining the series for the distinct 
families of solutions by an iterative process. 

As an example, we have applied this algorithm to obtain two term truncations 
of the series expansions of solutions to a cosmological model filled with a 
nonlinear causal viscous fluid. Thus, it is shown as feasible to deal with the 
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case branching that occurs in series solutions to nonlinear ordinary differential 
equations relevant to physical applications. 

It deserves to be investigated how the complexity of the algorithm increases 
with the order of iteration, and whether this growth puts an effective limit to 
practical calculations. In such it would be interesting to know whether 

more efficient algorithms could be devised. Also, it would be interesting to 
know whether expansions of solutions as shown in this paper can give infor- 
mation about the integrability of the equation. 
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